//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   FissionSource.cc
 * \author Jeremy Roberts
 * \date   Sep 4, 2011
 * \brief  FissionSource member definitions and explicit instantiations.
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//

#include "FissionSource.hh"

namespace slabtran
{

// Constructor.
template <class Equations, class D>
FissionSource<Equations, D>::FissionSource(SP_state state,
                                           SP_mesh mesh,
                                           SP_materials materials) :
  d_state(state),
  d_mesh(mesh),
  d_materials(materials),
  d_density(mesh, Equations::number_unknowns(), 0.0),
  d_source(mesh, Equations::number_unknowns(), 0.0),
  d_density_built(false),
  d_group(-1),    // negative useful for checks
  d_factor(-1.0)  // negative useful for checks
{
  //
}

// Fission source density.
template <class Equations, class D>
void FissionSource<Equations, D>::build_fission_density(double keff)
{
  // Store 1 / ( keff * norm ), where norm is 4*pi or 2.
  d_factor = Norm<D>::norm / keff;

  // Reset source density to zero.
  std::fill(d_density.begin(), d_density.end(), 0.0);

  // Grab cell material id's.
  const typename Mesh<D>::Vec_Id &matid = d_mesh->get_cell_id(MATERIAL_ID);

  // Loop through groups.
  for (int g = 0; g < d_state->number_groups(); g++)
  {

    // Grab the flux moments for this group.
    const Moments_Field_t &phi_g = d_state->get_phi(g);

    // Loop through spatial cells.
    for (int cell = 0; cell < d_mesh->number_cells(); cell++)
    {
      // Grab the cross sections for this cell.
      const typename Materials::XS &xs = d_materials->get(matid[cell]);

      // Grab nu and SigmaF for this group
      double nusig = xs.nu(g) * xs.sigma_f(g);

      // Loop through spatial unknowns in the cell and add to density.
      for (int n = 0, N = Equations::number_unknowns(); n < N; n++)
        d_density(cell, n) += nusig * phi_g(cell, n, 0, 0);

    } // end cell loop

  } // end group loop

  d_density_built = true;
}

// Fission group source.
template <class Equations, class D>
void FissionSource<Equations, D>::build_fission_source(int g)
{
  Require(d_density_built);

  // Store the group.
  d_group = g;

  // Grab cell material id's.
  const typename Mesh<D>::Vec_Id &matid = d_mesh->get_cell_id(MATERIAL_ID);

  // Loop through spatial cells.
  for (int cell = 0; cell < d_mesh->number_cells(); cell++)
  {
    // Grab the cross sections for this cell.
    const typename Materials::XS &xs = d_materials->get(matid[cell]);

    // Set factor * chi for this group
    double factor_chi = d_factor * xs.chi(g);

    // Loop through spatial unknowns in the cell and add to source.
    for (int n = 0, N = Equations::number_unknowns(); n < N; n++)
      d_source(cell, n) = factor_chi * d_density(cell, n);

  } // end cell loop

}

// Instantiations.
template class FissionSource<DD_1D_Equations, _1D> ;
template class FissionSource<DD_2D_Equations, _2D> ;
template class FissionSource<DD_3D_Equations, _3D> ;


} // end namespace slabtran

//---------------------------------------------------------------------------//
//              end of FissionSource.cc
//---------------------------------------------------------------------------//
